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A self-consistent method for calculating electron transport through a molecular device is devel- 
oped. It is based on density functional theory electronic structure calculations under periodic bound- 
ary conditions and implemented in the framework of the non-equilibrium Green function approach. 
To avoid the substantial computational cost in finding the I-V characteristic of large systems, we 
also develop an approximate but much more efficient non-self-consistent method. Here the change 
in effective potential in the device region caused by a bias is approximated by the main features of 
the voltage drop. As applications, the I-V curves of a carbon chain and an aluminum chain sand- 
wiched between two aluminum electrodes are calculated - two systems in which the voltage drops 
very differently. By comparing to the self-consistent results, we show that this non-self-consistent 
approach works well and can give quantitatively good results. 



PACS numbers: 73.40.Cg, 72.10.-d, 85.65.+h 



I. INTRODUCTION 

In recent years, electron transport through molecules 
sandwiched between metallic electrodes has been attract- 
ing increasing attention both for fundamental reasons 
and because it may form the basis of a future molecular 
electronics technologyi 2 ^^ Experimentally it is diffi- 
cult to precisely manipulate or even measure the atomic 
structure of the molecule-electrode contacts. Therefore, 
neither the influence of atomic structure on transport 
through the devices nor a path to improved performance 
is clear. As a result, the ability to calculate the atomic 
and electronic structure as well as the transport proper- 
ties of electrode-molecule-electrode systems is important 
and useful in this field. 

Electron transport through nanoscalc molecular de- 
vices differs significantly from that through macro- 
scopic semiconductor heterostructures. In the latter, the 
effective-mass approximation is generally successful be- 
cause of the periodic lattice structure and large elec- 
tron wavelength. In contrast, in a molecular device a 
carrier electron will be scattered by only a few atoms 
whose particular arrangement, then, matters a great 
deal. Consequently the effective-mass approximation 
breaks down, and the electronic structure of the molec- 
ular device must be taken into account explicitly. For 
this purpose, methods based on density functional theory 
(DFT) are sufficiently accurate and efficient Conven- 
tional DFT methods, however, deal with either closed 
molecular systems (in quantum chemistry) or periodic 
solids (in solid state physics), neither of which is appli- 
cable to molecular transport. Thus one needs to develop 
a DFT approach suitable for a system which is open, 
infinite, non-periodic, and non-equilibrium (if the bias 
voltage is nonzero). 

One way to do this was suggested by Lang, et al&^iSilL 
By using the jellium model for the two metallic elec- 
trodes of an electrode-molecule-electrode system, they 



mapped the Kohn-Sham equation of the system into the 
Lippmann-Schwinger scattering equation and solved for 
the scattering states self-consistently. They then cal- 
culated the current by summing up the contributions 
from all the scattering states, following a Landauer-type 
approach^ In this way, both the conductance and I-V 
characteristics of the system can be investigated. The use 
of the jellium model for electrodes is convenient and sim- 
ple but limited: it cannot include the effects of different 
contact geometries and surface relaxation, for instance. 
It also cannot deal with directional bonding such as in 
semiconductors and transition metals. As a result, the 
molecule-electrode charge transfer, which is one of the 
key factors affecting transport, may not be quantitatively 
correct^ 

Another way to develop the desired DFT approach 
is to use the non-equilibrium Green function (NEGF) 
method^^ The required open and non-equilibrium con- 
ditions can be treated rigorously, at least in a formal 
sense. This method is also closely related to the Lan- 
daur approach 12 and has proven to be powerful for study- 
ing electron transport through nanoscale devices. There- 
fore, by combining the NEGF method with conventional 
DFT-based electronic structure methods used in quan- 
tum chemistry or solid-state physics, the coherent trans- 
port properties of an electrode-molecule-electrode sys- 
tem can be determined fully self-consistently from first- 
principles. A further advantage of the NEGF+DFT com- 
bination is that the atomic structure of the device re- 
gion and the metallic electrodes are treated explicitly on 
the same footing. As has been mentioned, the molecule- 
electrode interaction will induce charge transfer between 
them and atomic relaxation of their contact - both have 
a significant effect on electron transport. As a result, 
the division of the system into the molecule and the elec- 
trodes is not meaningful anymore, and some parts of the 
electrodes must be included into the device region to form 
an "extended molecule" . 
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Based on this combined NEGF+DFT method, there 
are several successful implementationsiSiiLiSiiiLSl for 
molecular conduction and extensive theoretical results in 
the recent literaturei 16 i 17 i 18 i 19 i 20 i 21 i 22 i 23 According to the 
way of treating the extended molecule, the semi-infinite 
leads, and their couplings in a lead-molecule-lead (LML) 
system, these implementations can be roughly divided 
into two categories. 

In one category people adopted a cluster geometry for 
all the subsystems of a LML system or for the extended 
molecule with the leads treated by a tight binding ap- 
proach (for example, Ref. [THHH). It is then convenient 
to employ well-established quantum chemistry code (like 
Gaussian or DMol) to do the electronic structure calcu- 
lation for the subsystem(s). However, there are potential 
problems with these treatments for strong molecule-lead 
couplings: in this case it is obviously necessary to in- 
cluded large parts of the leads into the extended molecule 
so that the strong molecule-lead interaction can be fully 
accommodated. To eliminate the artificially introduced 
surface effects an even larger system is needed, which 
is usually difficult to deal with by a quantum chemistry 
code. So in practice only several (or even only one) lead 
atoms are attached to the m olecule to form an extended 
molecule (for example, Ref. I23l24j) . In this case, signif- 
icant artificial surface effects are inevitable, the contact 
atomic relaxation cannot be included, and an accurate 
molecule-lead coupling is not available. In addition, there 
may be artificial scattering at the interface between the 
tight-binding part of the lead and the DFT part of the 
lead (included in the extended molecule). 

In the other category (Ref ll8llT^) . people adopted pe- 
riodic boundary conditions (PBC) (as in solid state 
physics) with large parts of the leads included in the 
extended molecule, so that the interaction between the 
molecule and its images will be screened off by the metal- 
lic lead in between. In this case all the potential prob- 
lems mentioned above will be absent and the whole 
LML system becomes nearly perfect in geometry and all 
the subsystems are treated exactly on the same footing. 
Two examples of successful implementations adopting 
the PBC are the TranSiesta package 18 and the MCD- 
CAL packaged On the other hand, a drawback is intro- 
duced by PBC: when a bias is applied, the Hartree po- 
tential must jump unphysically between unit cells. This 
has previously been addressed by having an independent 
solution of the Poisson equation. 

In this paper, we first develop a fully self-consistent 
NEGF+DFT method with PBC, which has small but 
important differences from the two previous implemen- 
tations. The advantage of our method is that it is simple 
while still rigorous: the non-equilibrium condition under 
a bias is fully included in the NEGF part, and, as a result, 
we do not need to make changes in the conventional elec- 
tronic structure part. So it is straightforward to combine 
with any electronic structure method that uses a localized 
basis set. More importantly, in this way the problem of 
the unphysical jumps in the Hartree potential is avoided. 



A shortcoming of the full self-consistent (SC) 
NEGF+DFT approach is the large computational effort 
involved, especially for large systems, large bias volt- 
ages, or cases where many bias voltages need to be calcu- 
lated as for I-V characteristics. As a result, a non-self- 
consistent (non-SC) method with much higher efficiency 
and useful accuracy is highly desirable. As a step to- 
ward this goal, we also construct an approximate but 
much more efficient non-SC method in which the change 
in self-consistent effective potential in the device region 
caused by a bias is approximated by the main features of 
the voltage drop. 

As an application of our approach, in this paper we 
do calculations by combining it with a very efficient elec- 
tronic structure package SIESTAS The I-V curves of 
two systems with different typical voltage drop behaviors 
- a carbon or aluminum chain sandwiched between two 
aluminum electrodes - are calculated. Our self-consistent 
results are in good agreement with those from previous 
calculations By comparison to the self-consistent re- 
sults, we examine the validity of the non-SC approach, 
showing that this approach works quite well and can give 
quantitatively nearly correct answers. 

The arrangement of this paper is as follows. In Sec- 
tion II we give briefly a description of our implementation 
of the NEGF+DFT method. Because the basic formal- 
ism of the NEGF+DFT is well established, 16 ! 17 ! 18 ! 19 we 
show only those formula useful for introducing the new 
features of our method. The present SC and non-SC ap- 
proaches are explained in Section III. Section IV starts 
with results for a carbon chain and an aluminum chain 
sandwiched between two Al(OOl) electrodes. Our results 
are compared with previous results, and we discuss the 
validity of the non-SC approach by comparison to the 
self-consistent results. In Section V we summarize and 
conclude. 



II. NEGF+DFT METHOD AND ITS 
IMPLEMENTATION 

A. Modeling of real physical systems 

Experimentally, a molecular device system consists of 
at least a molecular junction coupled with two metallic 
electrodes (leads L and R) under a bias V\> (two-terminal 
system). In some cases, there is also a gate terminal ap- 
plying a gate voltage on the whole system (three-terminal 
system). Here we consider only the two-terminal system 
which is schematically shown in Fig. 1. An important 
consideration for modeling the real physical system is the 
charge transfer and atomic relaxation around the two 
molecule-lead contact regions. As a result, we have to 
include some parts of the metallic leads into the device 
region, forming an extended molecule. One obvious con- 
vergence criteria for the size of the extended molecule 
is its charge neutrality. Then the charge transfer and 
the potential disturbance caused by the molecule can be 
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considered screened off outside the extended molecule re- 
gion. In order to obtain good convergence, we actually 
include a large part of each metallic lead into the C re- 
gion, so that the layers adjacent to L and R (i.e., Cl 
and Cr parts in Fig. 1) have bulk properties. The total 
Hamiltonian of the system is: 



H = Hll + Hqc + Hrr + Hlc + Hcr- 



(1) 



Note that, here the leads L and R interact only through 
the molecular junction, so their direct interaction term 
Hlr vanishes (this can always be satisfied by using a 
localized basis set). 



Localized basis set 



When H is expanded in a basis set, generally only 
the matrix of Hcc (denoted Hcc) is finite. However, 
consider a localized (but not necessarily orthogonal) ba- 
sis set, by which we mean that the overlap between any 
two basis functions, (j>^(r — Ri) and </v(r — R2), will be 
zero if they are separated far enough from each other: 
S M „ = (fi\iy) = if |Ri — R2I > certain cutoff distance. 
In this case, the region C interacts directly only with 
finite parts of L and R, and the non-zero part of the ma- 
trices Hlc and Hcr also become finite. Furthermore, 
we can divide the leads L and R into principal layers so 
that any principal layer interacts only with its two near- 
est neighbors (see Fig. 1). As a result, the matrices Hll 



h' LL h' LL h'u, molecule 



h u L1 hV h°. 



R 



h RR h RR h Rg 



h jo h (J 
RF) 11 RR 11 RR 11 RR 



FIG. 1: Schematic drawing of a system containing a molecule 
sandwiched between two metallic electrodes (leads L and R). 
The region C is formed by including some parts of L and R so 
that the Cc part (extended molecule) is charge neutral and 
the Cl and Cr parts have bulk properties. Because of the 
use of a localized basis set, the leads L and R can be divided 
into principal layers (denoted by numbers 0, 1, 2, ...). Cl', 
Cqi , and Cri denote the parts used in the present non-SC 
approach (see Section II E). Their interface is called X in the 
text. h°'s and h^s are the Hamiltonian matrices within and 
between the principal layers, respectively. 



and Hrr have the following block tridiagonal form: 

if i - j = 
Hi- j = 1 

if j ~ i = 1 
if \i-j\ >1, 



(H LL )j,- 



(hi L ; 
0, 



it 



(2) 



where h° L and h LL are the Hamiltonian matrices within 
and between the principal layers, respectively, and i, j 
are principal layer indexes as shown in Fig. 1. Because 
the Cl and Cr parts which interact directly with L and 
R, have bulk properties, the non-zero part of Hlc (H-Cr) 
is just h\ L (h RR ), as shown in Fig. 1. 

In the localized basis and after the partition shown in 
Fig. 1 , the matrix Green function G of the whole system, 
defined by (ES - H)G(E) = I, satisfies 
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ES 



LL 
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rt 



LC 
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ESlc — Hlc 
H ( 
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(3) 



'CR 

Gll Glc Glr 
Gcl Gcc Gcr 

Grl Grc Grr 



The most important part of G is Gcc 7 corresponding to 
the region C; from the above equation, 

Gcc(E) = {ES CC - [Hcc + V L (E) + Xr(E)}}- 1 , 

(4) 

where Hl(E) and 'Er(E) are self-energies which incor- 
porate the effect of the two semi-infinite leads L and R, 
respectively. £l(_E), for example, is defined by 

H L (E) = (ESlc - H LC ) f G LL (E) (ES LC - Hlc) (5) 

where G LL is the retarded Green function of the left 
semi-infinite lead. The latter is given in turn by 

G° LL (E) = {zSll-HllY\ (6) 
z = E + ir], 

where a typical value for the lifetime broadening rj is 
about 1 meV. 

Because of the localized basis set, the non-zero part 
of S L c, H lc , &cr, and H C r become finite (being s\ L , 
h LL , s RR , and h RR ). As a result, only the part of G° LL 
and Gj^ corresponding to the th principal layer of the 
two leads (denoted g° LL and g° RR ) are needed for calcu- 
lating the non-zero part of the self-energies: 



V L (E) = (Esl L - hl L y g ° LL (E) (Es\ L hi L ) . (7) 

Our notation here follows that for the Hamiltonian: s 
and g are submatrices of the corresponding upper case 
matrices. g LL and g RR are simply the surface Green 
functions of the two semi-infinite leads. g° L , for example, 
can be calculated either by simple block recursion, 



sIl(E) = [zs 



LL 



n LL 



(8) 



(zsIl - ^ll) S°ll( e ) (zB 1 LL -h 1 LL )\ 
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or by a renormalization method 27 in terms of s^ L , s\ L , 
h° L , and \i\ L which can be determined by separate DFT 
calculations for the two leads. For small lifetime broaden- 
ing i] (1 meV), we find that the renormalization method 
is much faster than simple block recursion. This is nat- 
ural since n renormalization interations incorporate 2™ 
principal layers, while n recursions incorporate only n. 

From Gcc(E), the projected density of states (PDOS) 
on the molecule (indicated by m) is given by 

N m {E) = - -Im {Tr m [G cc {E + ir,) ■ S cc ] } , (9) 

where Tr m means the trace is performed only on the 
molecular part of the matrix. 



C. Current 

The Non-Equilibrium Green Function technique 
(NEGF)A£ii5i 2 &2& provides a convenient way to calculate 
the current by post-processing a DFT calculation. The 
result is quite natural and intuitive: First, the basic as- 
sumption is that there is no energy relaxation within the 
molecular region. Then, following a Landauer-like point 
of viewji 2 ^ one divides the electrons in the molecule into 
two sets using scattering-wave states, those that came 
from the left lead and those that came from the right. 
The left-lead states are, of course, filled up to the chem- 
ical potential in the left lead, p^, while the right-lead 
states are filled up to p R . In equilibrium, the two chem- 
ical potentials are equal, and the current carried by the 
left-lead states is, of course, equal to that carried by the 



right-lead states. As a bias is applied, the balance be- 
tween the two types of states is disrupted and current 
flows. As different states are populated because of the 
change in chemical potentials, the charge density in the 
molecule also changes. The potential profile must be 
solved for self-consistently in order to get an accurate 
measure of the transmission. It is this self-consistency 
which is the time-consuming part of the calculation. 

The expression for the steady-state current through 
the C region for applied bias Vb is 



2^ 
h 



/+oo 
T(E, V b ) [f(E-p L ) - f(E-p R )\ dE, 
-OO 



(10) 

where ph and p R are the chemical potentials, / is the 
Fermi function, and T(E, V&) is the transmission proba- 
bility for electrons from the left lead to right lead with 
energy E under bias Vb- The transmission probability is 
related to Green functions by 

T(E, V b ) = Tr [t l (E)Gcc(E)T r (E)G cc (E)] , (11) 

where 

T L . R {E) = i (X LM (E) - [H L , R (E)]^ (12) 
reflect the coupling at energy E between the C region 
and the leads L and R, respectively. 

The charge density corresponding to the above picture 
of left-lead states filled to p L and right-lead states filled 
p R can also be expressed in terms of Green functions. In 
particular, the density matrix of region C in the basis- 
function space is 



D 



cc 



1 r+°° r 

— / dE G C c(E)T L (E)G cc (E)f(E-p L )+G C c(E)T R (E)G cc (E)f(E-p R ) 

-- f + dElm[G CC (E)f(E-p L )] 

G cc (E)T R (E)Gl c (E)\ x [f(E - m ) - f(E - M£ )] . 

I 



(13) 
(14) 



2-J dE 



Time-reversal symmetry {G cc = G* cc ) was invoked in 
going from to l|14|) . The integrand of the first term of 
(|14fl is analytic (all poles of Gcc(E) are on real axis), so 
the integral can be evaluated easily by complex contour 
integration. However, the integrand of the second term is 
not analytic, so it must be evaluated by integrating very 
close to the real axis using a very fine energy mesh. The 
whole integration pathi 7 - is shown in Fig. 2. Because 
we construct the region C such that Cl and C R have 
essentially bulk properties, we can use the bulk density 
matrix for them. 



The calculated density matrix is then output to the 
DFT part to calculate the electron density p and to con- 
struct a new tlcc'- 

P (v) = [(D cc ) J ^(r), (15) 

(Hcc)m* = Mf + ^xt(r) + Mp(r)] + V«[p(r)] 



where T is the kinetic energy, and V cx t, Vh, and V^ c are 
the external, Hartree, and exchange-correlation potential 
energies, respectively. The new ~H.cc replaces the old, a 
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new Dec is calculated, and so on until Hcc or Dec 
converges. Finally, the transmission function T(E) can 
be calculated by Eq. ifTTl) . 

One subtlety here is the different boundary conditions 
used in the Green function and DFT parts - open versus 
periodic, respectively. This means that some iteration 
must be done even at V b = 0. If the supercell of the 
DFT part has the same size as region C in Fig. 1, then 
the Cl and Cr parts will interact directly due to the 
periodic boundary condition. However, this interaction 
is absent in the calculation of the density matrix. The 
same problem exists also for Hcc- So we need to do some 
translation work between the Green function and DFT 
parts: to add this interaction when we go from NEGF 
to DFT by using the density matrix elements between 
two adjacent principal layers, and to remove it when we 
go from DFT to NEGF by setting corresponding parts 
of Sec and Hcc to zero. Generally, the supercell of 
the DFT part can be made larger than the size of the 
region C, especially for systems without a translational 
symmetry, because the DFT part is usually much cheaper 
than the Green function part. 



III. NEW SELF-CONSISTENT AND 
NON-SELF-CONSISTENT APPROACHES 

A. Bias voltage 

For non-zero V b , care must be taken to account for the 
effects of the bias voltage on the charge density. One 
way to proceed is to apply a constant field in the direc- 
tion parallel to the leads within the supercell of the DFT 
calculation. Thus a linear drop is added to the external 
potential in Eq. i|16|) , and the effect on the charge density 
follows from, for instance, solving the Poisson equation. 
This approach is not straightforward for periodic bound- 
ary conditions because of the artificial potential jumps at 



Im[z] 



Re[z] 



Ml Mr 



FIG. 2: Schematic drawing of the integration path in the 
complex energy plane used to calculate the density matrix 
[Eq. C3J]. Eb is the lowest energy of occupied states, and 
fiL,R are the chemical potentials of the left and right leads, 
respectively < Hr is assumed). Note that for energy 
window [Eb,^l\ a complex contour integration is performed 
while for energy window \^ll,^r] a direct energy integration 
is performed by using a fine energy mesh and a very small 
imaginary part. 

the two supercell boundaries in the lead direction. One 
way to eliminate the unphysical jumps is to use a larger 
supercell for the DFT calculation and replace the Hamil- 
tonian of each of the regions near the potential jumps 
by the bulk one with a constant potential shift given by 
the bias voltage; this is implemented in the Transiesta 
program^ 

Here we propose a different approach to handle the 
bias, one which is less obvious but turns out to be sim- 
pler in the end: The bias is included through the density 
matrix (Dee) m the Green function calculation instead 
of the potential (Hec) m the DFT part. Specifically, 
we calculate the density matrix by Eq. I|13() under the 
boundary condition that there is a potential difference V b 
between part Cl (together with the left lead) and part 
Cr (together with the right lead). This is done by shift- 
ing all the potentials related to the left (right) lead and 
the C L (Cr) part by -V b /2 (+V b /2). Shifting the poten- 
tial in a lead is equivalent to directly shifting the energy 
by the opposite amount, so 



Z7T 



dE 



Gcc{E)T L {E + ^)G cc (£)/(£ - v L ) + G C c(E)Tr(E - e -^)G cc {E)f{E - m ) 



(17) 



Here the Green function Gcc{E) nas au the potential shifts included, 



Gcc(#) = ES 



H 



cc 



AH Ct +AH c „ + S L (£- 



*n(E-f) 



(18) 



r 



where the Cl and Cr parts of Hec are replaced by h^ L The actual computational process is exactly the same as 

that of Ijl4|l and Fig. 2, and self-consistency is achieved in 
the same way as for zero bias. Finally, in the calculation 



and h RR , and their potential shifts are 

[AH C J M „ = -y[Sd,„,V M ,,eC L , (19) 



of the transmission, the potential shift present in the self- 



fAH 



c R \ 



„ = +^[SccU,V/i,i/eCfl. (20) 
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energies appears in as in l|17|) . 

T(E,V b ) = T r {T L (E+^)G C c(E)T R (E-^)Gl c (E)} 

(21) 

While the two approaches mentioned here - the lin- 
ear external potential and the potential shift in the 
leads - give quite different results in the unphysical non- 
interacting limit, self-consistency ensures that they give 
the same result in physical cases. Our approach has the 
distinct technical advantage that the Green function and 
DFT calculations are completely disconnected, allowing 
the transport module to be easily combined with a wide 
variety of electronic structure codes. 



B. Approximate non-self-consistent approaches 

For large systems under large bias, the full SC ap- 
proach described above is computationally very difficult. 
The longest part of a one-time calculation is finding the 
surface Green functions (for all the points in the energy 
mesh), even with the fast renormalization method. How- 
ever, one only needs to do these calculations once and 
save the results. The major cost for a full self-consistent 
calculation is from finding G C c{E) by Eq. (@J at the 
many energies needed for the density matrix, especially 
for the very fine mesh in the energy window fi R ] (sec 
Fig. 2). 

To avoid the complex procedure and large computa- 
tional effort of full self-consistency, here we propose an 
approximate non-SC approach in which the bias is in- 
cluded by (a) applying a potential shift — eVb/2 (+eT4/2) 
to the left (right) lead through the self-energies as in 
Eq. I|21l) . and (b) introducing potential shifts non-self- 
consistently into the region C . The main consideration 
is that it may be a good approximation to replace the 
change in self-consistent effective potential caused by a 
bias by the main features of the voltage drop. This as- 
sumes, of course, that one can guess or motivate the main 
features in advance. For instance, for conductive molec- 
ular devices, the bias voltage will drop mainly at the left 
and right contact regions if the contacts have low trans- 
parency or over the whole molecular region if the contacts 
are very transparent. 

In our approach, we introduce new parts on the left 
(Cl 1 ) and right {Cr/) within the region C which extend 
from their respective leads to the molecular contacts, as 
shown in Fig. 1. We will denote by X the interface 
between the molecule Cc> and the regions Cl' plus Cri ■ 
If the voltage drop around the left (right) contact is aV b 
((1 — a)Vb), then the potential shifts are applied in the 
following way (which we call the AH1 treatment): 



Hcc + (° - 5) eV b s cc, (22) 
-eaV b [S C c}^,Vor veC L >, (23) 
[AH c ^] f = e{l-a)V h [S C c]^,liorv£C R ,(2±) 



Because the potential shift is applied to a matrix element 
when either orbital index is in the C'l',r' part, the voltage 
drop will be slightly smeared around the interface X. To 
explicitly show the role of Sec m l|23|l and 124fl . we also 
drop the potential using (called the AHO treatment) 



AH 



L -I {11/ 



-eaV b [Scc] M „ v G C L >, (25) 
[AHcJ^ = e(l-a)V b [S C cl„,»,v&C R ,. (26) 

In either case, Gcc{E) is determined by an equation 
analogous to Eq. I|18(l in the SC case, 



Gcc{E) — (^ES C c - 
+~Z L {E- 



H 



cc 



eVb, 
2 ' 



AH C - + AH 

eV b 



-V R (E 



) 



C, (27) 

-1 



The initial Hcc matrix here comes from a separate DFT 
calculation using a large L-C-R supercell. Finally, the 
current is calculated as usual through Eqs. ()21|l and 1)10(1 . 

So far we have considered the simplest case of low 
transparency contacts so that the voltage drops in only 
two places near the contacts. In this simplest case the 
a parameter here has the same role as the 77 parameter 
111 Ref. El For a real device, the voltage drop may be 
much more complicated, and there may be several dif- 
ferent voltage drops inside the device region. However, 
in our method all the factors affecting the voltage drop 
have been taken into account at the DFT level, and it 
is straightforward to generalize the present non-SC ap- 
proach for these more complicated situations (see calcu- 
lations for system B in Section IV) as long as the main 
features of the voltage drop are known. Actually, we will 
show later that results of this non-SC approach are not 
very sensitive to the choice of the voltage drop. If we 
assume that the form of the drop is not affected signif- 
icantly by a change in the bias voltage itself, then the 
main features of the voltage drop in a system can be de- 
termined by a single self-consistent calculation using a 
relatively small bias voltage. 



IV. APPLICATIONS: CHAINS OF CARBON OR 
ALUMINUM 

We report calculations of I-V curves for two systems: 
a carbon chain (system A) or an aluminum chain (system 
B) sandwiched between two aluminum leads in the (001) 
direction of bulk Al. The structures are shown in Fig. 3. 
No further atomic relaxation is performed for simplicity 
and for direct comparison with previous results. 

Our implementation of the transmission calculation is 
independent of the DFT part. Therefore, it can be eas- 
ily combined with any DFT package that uses a local- 
ized basis set. As an application, here we combine it 
with the very efficient full DFT package SIESTA,^ which 
adopts a LCAO-like and finite-range numerical basis set 
and makes use of pseudopotentials for atomic cores. In 
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FIG. 3: (Color online) Systems calculated: (a) a chain of 7 
carbon atoms sandwiched between two Al leads in the (001) 
direction of bulk Al, (b) a chain of 7 Al atoms sandwiched 
between the same leads. The C-Al distance in (a) is lA and 
the C-C distance is 1.323 A. The Al-Al (chain) distance in 
(b) is 2.86A and the Al-Al(surface) distance is 2.025A (i.e., 
the two Al atoms at the ends of the chain are in their bulk 
positions). The notations for different parts are the same 
as in Fig. 1. The numbers 0, 1, 2, 3, 4, 5, 6 denote different 
interfaces, called the interface X, between or C-r! an d Cqi 
which are used in the present non-self-consistent approach. 



our calculations we adopt a single zeta (SZ) basis set. 
To check the convergence of the results, we also calcu- 
late the equilibrium transmission function using a single 
zeta plus polarization (SZP) basis set. The difference is 
found to be minor. The PBE version of the generalized 
gradient approximation^ is adopted for the electron ex- 
change and correlation, and optimized Troullicr-Martins 
pseudopotentials^i are used for the atomic cores. The 
initial density matrix of the region C is obtained from a 
separate DFT calculation using a large L-C-R supercell. 

There are two main reasons for us to choose to study 
these systems. First, the transmission function T(E) of 
system A has been calculated by both the TranSiestai& 
and MCDCAL22 packages using a SZ basis set. So we 
can make a direct comparison to previous results. Sec- 
ond, systems A and B typify two different voltage drop 
behaviors (although both the C and Al chains are con- 
ductive). In system A, because the molecule-lead con- 
tact is a hetero-interface, the voltage drop will mainly 
occur around the two interfaces [see Fig. 4 (a), note that 
the voltage drops around the two interfaces are actually 
asymmetric]. In contrast, in system B the molecule-lead 
contact is a homo- interface, and furthermore the two Al 
atoms at the ends of the chain are at their bulk posi- 
tions. So the voltage drop will occur over the entire Al 
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FIG. 4: (Color online) Voltage drop in a plane going through 
the atomic chain for the Cc region of (a) the C chain and (b) 
the Al chain, as shown in Fig. 3, for an applied bias of -1 V. 
Note that (1) in (a) the voltage drop mainly occurs around 
the right contact region while in (b) it occurs over the entire 
chain, and (2) the oscillation in (b) is much larger than that 
in (a). This is because the electrons in the Al chain are more 
free than that in the C chain; as a result, the polarization 
induced by the bias is larger in the Al chain. 



chain in some way [see Fig. 4 (b)]. Our purpose is to see 
whether our non-SC approach can handle these different 
conditions. 



A. Transmission functions 

In Fig. 5 we show the calculated transmission func- 
tions and PDOS (projected on the chains) for system A 
and B under zero bias voltage. As it can be seen, the 
transmission function generally follows the PDOS except 
for some localized states (for instance, around 4 ~ 5 eV 
in (a) and 1 eV in (b)) which are not coupled with the 
left or the right lead. Our result of T(E) for system A 
is in very good agreement with the previous results from 
TranSiesta and MCDCAL packages!^ (see Fig. 6 (a) 
of Ref. E3). 

In Fig. 6 we show, for system A under a bias of 1.0V 
(a = 0.5), the calculated T(E) by the SC approach and 
the non-SC approach (AH1 treatment) with different 
choices for the interface X between the Cjy or Cr> and 
Cc" parts. The first thing we note is that our SC result 
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FIG. 5: Calculated transmission function T(E) (solid line) 
and projected density of states (PDOS, projected on the 
chains, dashed line) by the self-consistent approach for (a) the 
carbon chain (system A) and (b) the aluminum chain (system 
B) under zero bias voltage. When T(E) is significantly dif- 
ferent from the PDOS, it means that a localized state exists 
at that energy. 



for IV bias is in very good agreement with the previ- 
ous self-consistent resultsiS*^ (see Fig. 6 of Ref. list) . 
When we change the interface X (denoted by the differ- 
ent numbers in Fig. 6) from deep in the leads to the 
contact regions (i.e., interface X = 1 — > 2 — > 3), the non- 
SC result varies significantly. However, as X moves into 
the contact regions (i.e., interface X — 3, or 4, or 5), the 
result becomes very close to the SC result and insensi- 
tive to the exact position of X. This result is just what 
we expect because in system A the bias voltage drops 
mainly around the hetero-interface contact regions. This 
can be regarded as an advantage of the present non-SC 
approaches result is not strongly dependent on the tech- 
nical choice. 

Similar calculations of T(E) for system B under Vb = 
1.0V (a — 0.5, AH1 treatment for non-SC) are shown in 
Fig. 7. Again, after the interface X is moved into the 
contact regions the different non-SC result become quite 
close. However, compared to the case of system A, the 
agreement with the SC result is not good, especially in 
the energy range around the (averaged) Fermi level: the 
transmission from the non-SC calculations (except for in- 
terface X — 5) is noticeably larger than that from the SC 
calculation. This substantial disagreement originates in 
the difference between the self-consistent effective poten- 
tial [Fig. 4(b)] and the non-self-consistent one assumed 
in the non-SC approach. 
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FIG. 6: Transmission function T(E) for the carbon chain (sys- 
tem A) under 1 V bias (a = 0.5) from the SC approach and 
the non-SC approach (AH1 treatment) with different choices 
for the interface X between Cc> and the Cjy or Cri parts 
(denoted by '1', '5') as shown in Fig. 3 (a). Note the 
similarity between the different non-SC results and the fully 
SC result once X is placed in the contact regions (> 3). 



B. I-V curve of system A 

In order to show effects of different voltage drops and 
the difference between the AHO and AH1 treatments [see 
Eqs. (|23[l - i|2rj[) ] for potential shifting, we give in Fig. 8 
the calculated I-V curves for system A from the non-SC 
approach with interface X = 3 compared to the SC result. 
We do the non-SC calculations for three different voltage 
drops (a = 0.2, 0.5, 0.8) for the AH1 treatment. In ad- 
dition, for a = 0.8 we do the non-SC calculation with 
the AHO treatment. Among the three different voltage 
drops, the result for a = 0.2 is in poor agreement with 
the SC result while those for a = 0.5 and 0.8 are in 
good agreement. This makes sense in view of the main 
features of the voltage drop in Fig. 4(a): the bias volt- 
age will mainly drop around the left (right) contact for 
a positive (negative) bias. However, the small difference 
in I-V curve between a = 0.5 and 0.8 indicates that the 
I-V characteristics is actually not sensitive to the exact 
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FIG. 7: Calculated transmission function T(E) for the Al 
chain system under 1 V bias (a = 0.5, AH1 treatment for 
non-SC) by the SC and non-SC approaches. The notations 
are similar to those in Fig. 6 



change in voltage drop. By comparing the SC result and 
the two non-SC results for a — 0.8, we see that the AH1 
treatment improves the result markedly. This is under- 
standable because voltage drops in real physical systems 
are not sharp step-functions but are somewhat smeared 
out. 



C. I-V curves of systems A and B: comparison of 
different approaches 

I-V curves for systems A and B calculated by different 
approaches are given in Figs. 9 and 10, respectively. For 
the non-SC approach, a = 0.5 and the AH1 treatment 
are always adopted. It turns out that for small bias volt- 
ages (Vb < 0.3V) all the different treatments (the SC, 
the non-SC with different interface A), all give similar 
results, i.e., the effect of the different voltage drops is 
very small. Along with the increase of bias voltage, the 
difference among the different calculations becomes more 
and more significant. 

For the carbon chain (system A) , because the molecule- 
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FIG. 8: (Color online) Calculated I-V curves for the C chain 
system by the SC approach and the non-SC approach with 
interface X — 3 [as shown in Fig. 3 (a)]. Three different 
voltage drops are considered: a — 0.2, 0.5, 0.8. The AHO 
and AH1 treatments are adopted for a — 0.8. The simple 
non-SC approach does very well using a — 0.5 or 0.8 and the 
AH1 treatment. 



electrode contact is a hetero-interface and therefore the 




-1 1 
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FIG. 9: (Color online) Calculated I-V curves for the C chain 
system by the SC approach and the non-SC approach with 
different interface X indicated by the numbers, a = 0.5 and 
AH1 treatment are adopted for the non-SC calculations. 
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FIG. 10: (Color online) Calculated I-V curves for the Al chain 
system by the SC and non-SC approaches. The notations are 
similar to those in Fig. 9. '2+3+4+5+6' means a combined 
interface X (see Fig. 3 (b) ) and each interface bears 1/10 of 
a bias voltage (see Section III C). 

voltage drop occurs mainly around the contact regions, 
the SC result and the non-SC results with the interface 
X located around the two contact regions are in good 
agreement. As for the transmission function in Fig. 6, 
once the interface X is within the contact regions the re- 
sult is quite insensitive to the technical choice. This indi- 
cates that for systems made of conductive molecular junc- 
tions coupled with metallic electrodes through hetero- 
interfaces, the present non-SC approach works quite well 
and can give nearly quantitatively correct answers. 

For the aluminum chain (system B) the molecule- 
electrode contact is a homo-interface and the two Al 
atoms at the ends of the chain are at their bulk posi- 
tions; therefore, the voltage drop is not localized around 
the contact regions. Consequently, the result from the 
non-SC approach with the interface X located around 
the contact region is not in good agreement with the SC 
result. In order to further verify our analysis for sys- 
tem B, we generalize the present non-SC approach for a 
voltage drip occurring at multiple points: We use a com- 
bined interface X = 2+3+4+5+6 in which each layer 
bears a voltage drop of 14/10. Because of the role played 
by the overlap matrix in Eqs. (|23|l and Ij24(l . the result- 
ing voltage drop will occur over the entire device region. 
The calculated I-V curve by the generalized non-SC ap- 
proach is given in Fig. 10 by a violet solid line. The 
overall agreement with the SC result is remarkably im- 
proved, indicating that our analysis is reasonable. 



D. Limitation of the present approaches 

For finishing the discussion we would like to point out 
the cases where the present method will not work. Obvi- 
ously, the present method is only valid for steady state co- 
herent electron transport through metal-molecule-metal 
systems; there are basically two cases where our method 
does not work: (1) electron transport in Coulomb block- 
ade regime, for both the SC and non-SC approaches, and 
(2) cases where the main feature of voltage drop is sensi- 
tive to the value of the bias voltage itself, for the non-SC 
approach. In the first case, the contact barrier is so high 
that the molecule and the leads are essentially separated, 
and as a result, the molecular chemical potential is gen- 
erally different from the Fermi energies of the leads even 
under zero bias. Because in our DFT+NEGF approach 
there is only one Fermi energy under zero bias, it will 
fail in this case. The second case is just the opposite to 
that assumed in our non-SC approach. We don't know 
at this moment what systems will have this behavior or 
whether such kind of systems exist. But this can be eas- 
ily checked by doing selfconsistent calculations for several 
different bias voltages within the bias range interested. 



V. SUMMARY 

A full self-consistent DFT-electronic-structure-based 
Green function method has been proposed and imple- 
mented for electron transport from molecular devices. 
Our method is simple and straightforward while strict. 
The implementation is very independent of the DFT elec- 
tronic structure part; it can be easily combined with 
any electronic structure package using a localized ba- 
sis set. In an effort to avoid the extremely burdensome 
computational cost for large systems or for I-V charac- 
teristic analysis, we developed an approximate non-self- 
consistent approach in which the change in effective po- 
tential caused by a bias in the device region of a system is 
approximated by the main features of the voltage drop. 

As applications of our methods, we calculated the I- 
V curves for two different systems with different typical 
voltage drops: a carbon chain and a aluminum chain 
sandwiched between two aluminum electrodes. Our self- 
consistent results are in very good agreement with those 
from other calculations. For both systems the present 
non-SC approach can give results in good agreement with 
the self-consistent results, indicating that it is a good ap- 
proximate method with high efficiency for I-V character- 
istic analysis^ (more then one order of magnitude faster 
for moderate systems). It is straightforward to general- 
ize this non-SC approach to deal with any kind of voltage 
drop situation. 

This work was supported in part by the NSF (DMR- 
0103003). 
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